#! /usr/bin/env python

'''
This script simulates a neuron driven by an excitatory and an
inhibitory population of neurons firing Poisson spike trains. The aim
is to find a firing rate for the inhibitory population that will make
the neuron fire at the same rate as the excitatory population.

Optimization is performed using the bisection method from Scipy,
simulating the network repeatedly.

This example is also shown in the article Eppler et al. (2009) PyNEST:
A convenient interface to the {NEST} simulator, Front. Neuroinform.
doi:10.3389/neuro.11.012.2008
'''

# First, we import all necessary modules for simulation, analysis and
# plotting. Additionally, we set the verbosity to suppress info
# messages

from scipy.optimize import bisect

import nest
import nest.voltage_trace

nest.sr("M_WARNING setverbosity")
nest.ResetKernel()

# Second, the simulation parameters are assigned to variables.

t_sim = 25000.0    # how long we simulate
n_ex  = 16000      # size of the excitatory population
n_in  =  4000      # size of the inhibitory population
r_ex  =     5.0    # mean rate of the excitatory population
r_in  =    20.5    # initial rate of the inhibitory population
epsc  =    45.0    # peak amplitude of excitatory synaptic currents
ipsc  =   -45.0    # peak amplitude of inhibitory synaptic currents
d     =     1.0    # synaptic delay
lower =    15.0    # lower bound of the search interval
upper =    25.0    # upper bound of the search interval
prec  =     0.01   # how close need the excitatory rates be

# Third, the nodes are created using Create(). We store the returned
# handles in variables for later reference.

neuron        = nest.Create("iaf_neuron")
noise         = nest.Create("poisson_generator",2)
voltmeter     = nest.Create("voltmeter")
spikedetector = nest.Create("spike_detector")

# Fourth, the excitatory Poisson generator (noise[0]) and the
# voltmeter are configured using SetStatus, which expects a list of
# node handles and a list of parameter dictionaries. The rate of the
# inhibitory Poisson generator is set later. Note that we need not set
# parameters for the neuron and the spike detector, since they have
# satisfactory defaults.

nest.SetStatus(noise, [{"rate": n_ex*r_ex}, {"rate": n_in*r_in}])
nest.SetStatus(voltmeter, {"interval": 10.0, "withgid": True, "withtime": True})

# Fifth, the neuron is connected to the spike detector and the
# voltmeter, as are the two Poisson generators to the neuron. The
# command Connect() has different variants. Plain Connect() just takes
# the handles of pre- and post-synaptic nodes and uses the default
# values for weight and delay. ConvergentConnect() takes four
# arguments: A list of pre-synaptic nodes, a list of post-synaptic
# nodes, and lists of weights and delays. Note that the connection
# direction for the voltmeter is reversed compared to the spike
# detector, because it observes the neuron instead of receiving events
# from it. Thus, Connect() reflects the direction of signal flow in
# the simulation kernel rather than the physical process of inserting
# an electrode into the neuron. The latter semantics is presently not
# available in NEST.

nest.Connect(neuron, spikedetector)
nest.Connect(voltmeter, neuron)
nest.ConvergentConnect(noise, neuron, [epsc, ipsc], 1.0)

# To determine the optimal rate of the neurons in the inhibitory
# population, the network is simulated several times for different
# values of the inhibitory rate while measuring the rate of the target
# neuron. This is done until the rate of the target neuron matches the
# rate of the neurons in the excitatory population with a certain
# accuracy. The algorithm is implemented in two steps:
# 
# First, the function output_rate() is defined to measure the firing
# rate of the target neuron for a given rate of the inhibitory
# neurons.

def output_rate(guess):
    print "Inhibitory rate estimate: %5.2f Hz ->" % guess,   
    rate = float(abs(n_in*guess))
    nest.SetStatus([noise[1]], "rate", rate)
    nest.SetStatus(spikedetector, "n_events", 0)
    nest.Simulate(t_sim)
    out=nest.GetStatus(spikedetector, "n_events")[0]*1000.0/t_sim
    print "Neuron rate: %6.2f Hz (goal: %4.2f Hz)" % (out, r_ex)
    return out

# The function takes the firing rate of the inhibitory neurons as an
# argument. It scales the rate with the size of the inhibitory
# population and configures the inhibitory Poisson generator
# (noise[1]) accordingly. Then, the spike-counter of the spike
# detector is reset to zero. The network is simulated using
# Simulate(), which takes the desired simulation time in milliseconds
# and advances the network state by this amount of time. During
# simulation, the spike detector counts the spikes of the target
# neuron and the total number is read out at the end of the simulation
# period. The return value of output\_rate() is the firing rate of the
# target neuron in Hz.
# 
# Second, the scipy function bisect() is used to determine the optimal
# firing rate of the neurons of the inhibitory population.

in_rate = bisect(lambda x: output_rate(x)-r_ex, lower, upper, xtol=prec)
print "\nOptimal rate for the inhibitory population: %.2f Hz" % in_rate

# The function bisect() takes four arguments: First a function whose
# zero crossing is to be determined. Here, the firing rate of the
# target neuron should equal the firing rate of the neurons of the
# excitatory population. Thus we define an anonymous function (using
# lambda) that returns the difference between the actual rate of the
# target neuron and the rate of the excitatory Poisson generator,
# given a rate for the inhibitory neurons. The next two arguments are
# the lower and upper bound of the interval in which to search for the
# zero crossing. The fourth argument of bisect() is the desired
# relative precision of the zero crossing.
# 
# Finally, we plot the target neuron's membrane potential as a
# function of time.

nest.voltage_trace.from_device(voltmeter)
nest.voltage_trace.show()
